Kristin Broms, Neptune & Co.
October 18, 2017
Static maps
3D maps
Interactive maps
Animated maps
Areal data exploration
Additional packages
As you can see, “spatial data” has come to mean much more than creating a map! Visualizing spatial data is an active area of software development, with many exciting packages and functions becoming available on a regular basis.
library(rgdal)
# streams <- readOGR(dsn = "../data/R_GIS_Layers/",
# layer = "Cove_Drainage_WGS84")
summary(streams)## Object of class SpatialLinesDataFrame
## Coordinates:
## min max
## x -109.35383 -109.12624
## y 36.46573 36.94164
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## OBJECTID ComID GNIS_ID GNIS_Name
## Min. : 0 Min. : 0 00003435: 2 Cove Wash: 2
## 1st Qu.: 32992 1st Qu.:103664065 NA's :47 NA's :47
## Median : 65593 Median :103665843
## Mean : 68116 Mean :101439114
## 3rd Qu.:105309 3rd Qu.:126660550
## Max. :136421 Max. :126669825
##
## ReachCode FType FCode Source
## 14080204001554: 2 Min. :460 Min. : 0 DOQQ: 1
## 14080204001982: 2 1st Qu.:460 1st Qu.:46003 NHDH:48
## 14080105001404: 1 Median :460 Median :46003
## 14080105001406: 1 Mean :462 Mean :45264
## 14080105001455: 1 3rd Qu.:460 3rd Qu.:46003
## (Other) :41 Max. :558 Max. :55800
## NA's : 1
## Name
## Middle 1 : 2
## Cove Wash : 1
## Cove Wash Middle: 1
## Cove Wash North : 1
## Cove Wash South : 1
## (Other) :17
## NA's :26
plot(streams, col="blue")## note that plots using base R may be distorted.library(ggplot2)
streams_gg <- fortify(streams) ## or use broom::tidy()
# str(streams_gg)
ggplot() +
geom_path(data = streams_gg,
aes(x = long, y = lat, group = group),
color = "blue", size = 1.5)## note that without a coordinates correction, the plot may be distorted.library(ggmap)
# Might need older version of ggplot2:
# devtools::install_github("hadley/ggplot2@v2.2.0")
# library(ggplot2)
myLocation <- c(lon = -109.2279, lat = 36.545)
myMap <- get_map(location = myLocation, source = "google",
maptype = "satellite", crop = FALSE, zoom = 13)
g <- ggmap(myMap, darken = c(0.25, "white"))
g +
geom_path(data = streams_gg,
aes(x = long, y = lat, group = group),
color = "lightblue", size = 1.5)
ggsave("shp_and_ggmap.png", dpi = 72) library(dplyr)
alum_data <- plot_data %>% filter(Analyte == "ALUMINUM")
g +
geom_path(data = streams_gg,
aes(x = long, y = lat, group = group),
color = "lightblue", size = 1.5) +
geom_point(data=alum_data,
aes(x = Longitude, y = Latitude, size = Result, fill = Result),
shape = 21, alpha = 0.9) +
scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
scale_size(range = c(2, 8)) +
guides(size = FALSE) +
coord_equal()
ggsave("static_map.png", dpi = 72) head(alum_data)## Analyte Result Units Latitude Longitude
## 1 ALUMINUM 9120.859 mg/kg 36.59625 -109.1736
## 2 ALUMINUM 7429.408 mg/kg 36.59625 -109.1736
## 3 ALUMINUM 1657.721 mg/kg 36.61342 -109.1356
## 4 ALUMINUM 1709.072 mg/kg 36.61342 -109.1356
## 5 ALUMINUM 5867.922 mg/kg 36.56392 -109.2120
## 6 ALUMINUM 46050.939 mg/kg 36.56392 -109.2120
Common ArcGIS functions, such as unions and intersections, can be done in R. The main advantages are that all manipulations are recorded and reproducible, and code can be automated.
library(scatterplot3d)
with(alum_data, {
lollipop_plot <- scatterplot3d(streams_gg$long, streams_gg$lat,
rep(0, length(streams_gg$long)), # x y and z axis
color = "blue", pch = 16, type = "p",
angle = 120,
scale.y = 1.75,
main = "Example Lollipop plot",
xlab = "Longitude",
ylab = "Latitude",
zlab = "Concentration values",
xlim = c(-109.273, -109.15),
ylim = c(36.50, 36.59),
zlim = range(0, 46000),
grid = TRUE,
box = FALSE)
# add the legend
legend("topright", inset = 0.07,
bty = "n",
title = "Type of Results",
c("High","Low"),
fill = c("#1B9E77", "#D95F02"))
# add the lollipop points
lollipop_plot$points3d(Longitude, Latitude, Result,
col = "#282830",
pch = 21,
bg = ifelse(Result > 10000, "#1B9E77", "#D95F02"),
lwd = 1,
type = "h",
cex = (Result / 50000) + 1)
})lollipop plot
(Generated using a modified version of the scatterplot3d function: https://github.com/USEPA/R-User-Group/tree/master/contributedCode/scatterplot3dMap)
Reference: Beaulieu, et al. 2016. Estimates of reservoir methane emissions based on a spatially balanced probabilistic-survey. Limnology and Oceanography, 61: S27-S40.
library(rgl)
library(dplyr)
plot_3d <- with(alum_data,
plot3d(Longitude, Latitude, log(Result), # x y and z axis
col = ifelse(Result > 10000, "#1B9E77", "#D95F02"), size = 5,
type = "p"))
rglwidget(elementId = "plot3drgl") # to show in presentation## plotly recommends the development version of ggplot2, but will also work with the
## latest version of the package, version 2.2.1
# library(devtools)
devtools::install_github('hadley/ggplot2')
library(ggplot2)
library(plotly)
p <- plot_ly(alum_data,
x = ~Longitude, y = ~Latitude, z = ~Result,
marker = list(color = ~Result,
colorscale = c('#FFE1A1', '#683531'), showscale = TRUE)) %>%
add_markers() pinteractive_map <- g +
geom_point(data = alum_data,
aes(x = Longitude, y = Latitude, size = Result, fill = Result),
color = "#000000", shape = 21, alpha = 0.8) +
scale_fill_gradient(low = "#1f77b4", high = "#d62728",
guide = "legend") +
scale_size(range = c(2, 8)) +
coord_equal() +
ggtitle("ALUMINUM")ggplotly(interactive_map + theme_void(), filename="plotly")ggmap and plotly are not completely compatible. (ggmap uses ggplot2.2.0, plotly wants the latest development version of ggplot2.)
“ggiraph” also makes ggplot figures interactive. It currently has limited uses but may be expanding.
library(animation)
ani.record(reset = TRUE)
for(a in unique(plot_data$Analyte)) {
plot_data_a <- subset(plot_data, Analyte == a)
animated_maps <- g +
geom_path(data = streams_gg,
aes(x = long, y = lat, group = group),
size = 1.5, color = 'lightblue') +
geom_point(data = plot_data_a,
aes(x = Longitude, y=Latitude, size = Result, fill = Result),
shape = 21, alpha = 0.8) +
scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
scale_size(range = c(2, 8)) +
guides(size = FALSE) +
coord_equal() +
ggtitle(a)
print(animated_maps)
ani.record()
}
oopts = ani.options(interval = 1)
ani.replay()
saveHTML(ani.replay(), img.name = "animation_plot",
htmlfile = "animation_5metals.html", ani.width=600, ani.height=400)Alternate package: library(gganimate)
Alternate package: library(anim.plots)
# Show neighbors on map of subbasins
library(sp)
class(lulc.sp1)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
spplot(lulc.sp1, zcol = "PCTHERB", col = "black",
main = "Percent Herb")library(spdep)
## determine who shares a border
ctchmt.nb1 <- poly2nb(lulc.sp1, queen = TRUE)
# put neighbors into list
ctchmt.nbwts.list <- nb2listw(ctchmt.nb1, style = "B")
moran.test(lulc.sp1$PCTHERB, ctchmt.nbwts.list)##
## Moran I test under randomisation
##
## data: lulc.sp1$PCTHERB
## weights: ctchmt.nbwts.list
##
## Moran I statistic standard deviate = 2.9826, p-value = 0.001429
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.212182899 -0.016949153 0.005901866
par(mar=c(5, 5, 2, 2))
moran.plot(lulc.sp1$PCTHERB, listw = ctchmt.nbwts.list,
labels = lulc.df1$PCTHERB,
xlab = "PCTHERB", ylab = "Spatially Lagged PCTHERB",
main = "Moran Plot for PCT HERB")## The network graph may require the following graph package that need to be
## installed from source. (This package is not generally required to create
## similar graphs; they are only required for the example. See the tutorial for
## other examples based on csv files.)
# source("https://bioconductor.org/biocLite.R")
# biocLite("graph")
# biocLite("Rgraphviz")
library(magrittr) ## for %<>% function
library(visNetwork) ## to create the interactive network graph
library(jsonlite)
library(igraph)
network_interactive <- function(g, nodes, coords){
## first, switch graph construction to data frames:
nNodes <- length(nodes)
nodes_df <- matrix(NA, nrow=nNodes, 4) # save first 4 attributes
for (i in 1:nNodes){
nodes_df[i, ] <- unlist(jsonlite::fromJSON(nodes[[i]]))[1:4]
}
nodes_df <- as.data.frame(nodes_df, stringsAsFactors=FALSE)
## because edges uses the names, need id to equal names of nodes, not numbers
names(nodes_df)[1:3] <- c("old_id", "id", "type")
## and id column needs to go first: (only for igraph package, not visNetwork)
## igraph packge was used to help decide node coordinates
nodes_df <- nodes_df[, c('id', 'old_id', 'type')]
g_edges <- Rgraphviz::buildEdgeList(g)
nEdges <- length(g_edges)
edges_df <- NULL
for (i in 1:nEdges){
tmp <- c(from=g_edges[[i]]@from,
to=g_edges[[i]]@to,
unlist(g_edges[[i]]@attrs))
edges_df %<>% bind_rows(tmp)
}
edges_df$weight <- as.numeric(edges_df$weight)
## Add layer info/ color column
nodes_df$col_layer <-
ifelse(nodes_df$id %in% grep("MainInput", nodes_df$id, value=T), 1, #
ifelse(nodes_df$id %in% grep("OtherInput", nodes_df$id, value=T), 3,
ifelse(nodes_df$id %in% grep("Midvalue", nodes_df$id, value=T), 4,
ifelse(nodes_df$id %in% grep("MidEqn", nodes_df$id, value=T), 5,
6))))
## add coordinates to determine layout of each node:
nodes_df <- full_join(nodes_df, coords)
## specify colors to use for each col_layer
graph_colors <- c("gold", "darkorange", "tomato",
"palegreen", "seagreen", "royalblue")
# frame = borders
frame_colors <- c("gold3", "darkorange3", "tomato4",
"palegreen3", "seagreen4", "royalblue4")
## (visNetwork doesn't like color names with numbers)
frame_colors_rgb <- rgb(t(col2rgb(frame_colors)), maxColorValue = 255)
## Add attributes so that the graph looks good:
nodes_df$shape <- "ellipse"
nodes_df$color.background <- graph_colors[nodes_df$col_layer]
nodes_df$color.border <- frame_colors_rgb[nodes_df$col_layer]
nodes_df$color.highlight.border <- "darkred"
edges_df$arrows <- "to" # draw arrows on the edges.
edges_df$color.highlight <- "black"
network <- visNetwork::visNetwork(nodes_df, edges_df, width="100%", physics = TRUE) %>%
visNetwork::visIgraphLayout() %>%
visNetwork::visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE,
manipulation = FALSE) %>%
visNetwork::visEdges(color = "slategray", smooth = list(roundess = 1)) %>%
visNetwork::visNodes(font = list(size = 22), shape="ellipse")
return(network)
}network_interactive(g = networkDag, nodes = networkNodes, coords = networkCoords)par(mar=rep(0, 4)) ## change margins so plot fills entire space
## default colors:
plot(1:8, pch=16, cex=5, col=1:8)## Defining a new palette
new_palette <- c("darkred", "chartreuse", "turquoise", "purple",
"gray45", "plum", "black", "#F08800")
palette(new_palette)
plot(1:8, pch=16, cex=5, col=1:8)palette("default") ## return palette to default colorssessionInfo()## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] igraph_1.1.2 jsonlite_1.5 visNetwork_2.0.1
## [4] magrittr_1.5 spdep_0.6-15 Matrix_1.2-11
## [7] bindrcpp_0.2 plotly_4.7.1 dplyr_0.7.4
## [10] rgl_0.98.1 ggplot2_2.2.1.9000 rgdal_1.2-13
## [13] sp_1.2-5
##
## loaded via a namespace (and not attached):
## [1] httr_1.3.1 maps_3.2.0 tidyr_0.7.1
## [4] viridisLite_0.2.0 splines_3.4.2 gtools_3.5.0
## [7] shiny_1.0.5 assertthat_0.2.0 expm_0.999-2
## [10] stats4_3.4.2 yaml_2.1.14 LearnBayes_2.15
## [13] backports_1.1.1 lattice_0.20-35 glue_1.1.1
## [16] digest_0.6.12 colorspace_1.3-2 htmltools_0.3.6
## [19] httpuv_1.3.5 plyr_1.8.4 pkgconfig_2.0.1
## [22] devtools_1.13.3 gmodels_2.16.2 purrr_0.2.3
## [25] xtable_1.8-2 scales_0.5.0.9000 gdata_2.18.0
## [28] jpeg_0.1-8 ggmap_2.6.1 git2r_0.19.0
## [31] tibble_1.3.4 withr_2.0.0 BiocGenerics_0.22.1
## [34] lazyeval_0.2.0 proto_1.0.0 mime_0.5
## [37] deldir_0.1-14 memoise_1.1.0 evaluate_0.10.1
## [40] nlme_3.1-131 MASS_7.3-47 graph_1.54.0
## [43] tools_3.4.2 data.table_1.10.4-1 geosphere_1.5-5
## [46] RgoogleMaps_1.4.1 stringr_1.2.0 munsell_0.4.3
## [49] compiler_3.4.2 rlang_0.1.2 grid_3.4.2
## [52] rjson_0.2.15 htmlwidgets_0.9 crosstalk_1.0.0
## [55] base64enc_0.1-3 labeling_0.3 rmarkdown_1.6
## [58] boot_1.3-20 gtable_0.2.0 curl_3.0
## [61] reshape2_1.4.2 R6_2.2.2 knitr_1.17
## [64] bindr_0.1 rprojroot_1.2 Rgraphviz_2.20.0
## [67] stringi_1.1.5 parallel_3.4.2 Rcpp_0.12.13
## [70] mapproj_1.2-5 png_0.1-7 coda_0.19-1